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Preface 


IMs  study  began  with  an  attempt  to  explain,  theoretically,  or 
experimentally,  how  the  Increased  rate  of  convergence  was  achieved 
for  two  methods  of  iteration  previously  developed  at  the  Air  Force 
Institute  of  Technology.  Analysis  of  the  simpler  method  applied  to 
point  successive  overrelaxation  was  attempted  first.  Difficulties  in 
applying  the  theory  to  the  method  and  achieving  the  reported  increase 
In  convergence  rates  in  actual  test  runs,  were  immediately  encountered. 

The  length  of  time  needed  to  finally  resolve  these  difficulties  pre¬ 
cluded  analysis  of  the  other  method,  which  was  applied  to  block 
successive  overrelaxation. 

Approximately  five  hours  of  oomputer  time  on  an  IK  I  7094  com¬ 
puter  were  used  in  testing  and  comparing  various  scanning  techniques 
with  the  method  previously  developed.  The  method  of  point  successive 
overrelaxation  was  applied  to  the  five  point  finite  difference  approx¬ 
imations  to  Laplace's  equation  and  the  transient  heat  transfer 
equation. 

Two  scanning  methods,  which  I  have  elected  to  call  the  standard 
scan  and  the  even-odd  scan,  were  programed  along  with  the  method 
previously  developed.  It  was  found  that  the  standard  and  even -odd 
scans  were  iteratively  faster  as  the  number  of  required  iterations 
was  increased.  Both  of  these  soans  are  "consistent"  and  therefore 
qualify  for  the  maximum  asymptotic  rate  of  convergence  amongst  the 
group  of  all  possible  scans.  I  have  attempted  in  chapter  II  to  ex¬ 
plain  explicitly  how  one  determines  when  a  given  scan  is  consistent. 

The  method  Is  due  to  Varga  and  can  be  easily  applied  once  It  is  understood. 
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Abstract 

This  study  investigates  the  affect  of  ordering  of  finite  differ¬ 
ence  equations  on  the  average  rate  of  convergence  of  the  point  success¬ 
ive  overrelaxation  iterative  method.  The  method  is  applied  to  the  two 
dimensional  five  point  implicit  finite  difference  representation  of 
Laplace's  equation  and  the  diffusion  equation  of  transient  heat 
transfer. 

The  average  rate  of  convergence  for  a  specified  number  of  itera¬ 
tions  is  found  to  depend  on  the  magnitude  of  the  projection  of  the 
initial  error  vector  onto  the  dominant  eigenvector  of  the  point  suc¬ 
cessive  overrelaxation  iteration  matrix.  Although  the  asymptotic  rate 
of  convergence  has  been  proven  to  be  the  same,  and  optimum,  for  all 
consistent  orderings,  the  average  rate  of  oonvergenoe  for  a  specific 
number  of  iterations  is  shown  to  be  different  for  two  consistent 
orderings.  For  a  redaction  in  the  Initial  error  vector  by  three 
orders  of  magnitude,  a  difference  of  20#  in  average  convergence  rates 
is  achieved.  An  increase  in  the  orders  of  magnitude  reduction  pro¬ 
duces  a  proportionate  decrease  in  the  per  cent  difference  between  the 
two  consistent  scans. 
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I.  Introduction 

With  the  growth  in  speed  and  complexity  of  modern  high  speed 
digital  computers  has  come  an  Increase  in  interest  in  finding  or 
approximating  the  solution  to  partial  differental  equations  in  several 
variables  by  use  of  these  machines.  It  is  naturally  desirable  to  find 
a  method  which  will  yield  a  solution  to  the  matrix  problems  arising 
from  discrete  approximations  to  partial  differential  equations  in  the 
mimimura  time  on  the  computer.  This  will  allow  more  and/or  larger 
problems  to  be  solved  with  the  required  accuracy. 

From  the  early  work  of  Young  (Ref  7)  and  Frankel  (Ref  4)  in  1950, 
the  method  of  point  successive  overrelaxation  (to  use  Young's  nomenclature) 
has  been  expanded  to  cover  a  wide  range  of  matrix  equations  through  the 
concept  of  p-cyclic  matrices  treated  by  Varga  (Ref  6).  This  theory 
also  applies  to  the  newer  method  of  block  successive  overrelaxation. 

The  mathematical  theory,  although  it  does  provide  a  good  general  frame¬ 
work  and  basis  for  direction  in  research,  is  not  sufficiently  ooapleta 
to  determine  tne  optimum  method. 

For  large  problems  the  most  successful  technique  to  date  employ* 
the  use  of  an  implicit  approximation  to  the  differential  equation 
wherein  the  correfronr) ing  difference  equations  *r->  col  red  repeatedly 
until  two  successive  iterations  yield  solutions  which  differ  only  by 
a  small  predetermined  amount.  The  discrete  values  so  obtained  are 
then  taken  as  the  solution  of  the  differential  equation  at  the 
respective  points.  Various  methods  have  bean  developed  which  use 
this  general  approach. 
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Previous  work  at  the  Air  Force  Institute  of  Technology  by  Cudahy 
(1965)  (Ref  2)  with  point  successive  overrelaxation  and  Wright  (1967) 
(Ref  10)  with  block  successive  overrelaxation  showed,  by  an  essentially 
empirical  approach,  that  the  ordering  of  mes v  points,  or  blocks  of 
mesh  points,  for  solution  of  the  difference  equations  can  have  a 
significant  effect  on  the  rate  of  convergence  of  the  iterative  method. 

Purpose 

The  objective  of  the  present  study  was  to  investigate  the  effect 
of  ordering,  or  method  of  mesh  point  scanning,  on  the  average  and 
asymptotic  rates  of  convergence  of  the  point  successive  overrelaxation 
iterative  method,  and  to  account  for  the  differences  in  experimentally 
observed  rates  of  convergence  for  different  scanning  methods. 

Method 

The  point  successive  overrelaxation  ito  ive  method  was  applied 
to  the  five  point  finite  difference  approximations  to  Laplace's  equa¬ 
tion  and  the  heat  equation.  The  equations  were  applied  over  a  square 
region  in  two  dimensions  so  that  the  optimum  relaxation  factor  could 
be  pi adicted  theoretically  for  Laplace's  equation,  and  so  that  the  pre¬ 
sent  data  could  be  correlated  with  the  previous  work  of  Cudahy  (Ref  2) 
for  the  heat  equation.  The  ability  to  predict  the  optimum  relaxation 
factor  greatly  reduces  the  compif  r  time  necessary  tc  obtain  the  desired 
data.  For  the  heat  equation  solutions  were  obtained  with  several 
relaxation  factors  for  each  t- ><><■<  step.  The  optimum  relaxation  factor 
is  different  for  each  time  step. 

The  even-odd  scan  (denoted  odd-even  parity  by  Forsythe  and  Wasow) 
(Ref  3:to0)  and  the  method  developed  by  Cudahy  (Ref  2:3*0  were  compared 
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with  the  standard  scan  used  by  Young  (Ref  8).  The  results  of  those 
comparison^  are  presented  In  chapter  Ill.  Chapter  IV  contains  a  dis¬ 
cussion  of  the  results  of  this  study  in  terms  of  the  theory  presented 
In  chapter  II. 
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II.  Theory 


Point  Successive  Overrelaxatlon 

To  derive  the  point  successive  overralaxation  iterative  method 
consider  a  system  o f  linear  equations 


ki» 


1  *  i  *  n 


(1) 


where  the  a^  j  are  elements  of  the  n  x  n  complex:  matrix  A.  This  can 
be  written  in  matrix  notation  as 

A  x  =  k  (2) 

where  k  is  a  given  column  vector.  The  solution  to  this  equation  exists 
and  is  unique  if  and  only  if  the  matrix  A  is  nonsingular  (Ref  6:56). 

The  solution  vector  can  then  be  written  explicitly  as 

A”1  k  (3) 


The  diagonal  elements  *±9±  of  A  are  now  assumed  to  be  nonzero  complex 
numbers - 

The  matrix  A  can  bf  'scprossed  os 

A  «  D  -  S  -  ?  (4) 


where  D  is  a  strictly  diagonal  matrix,  E  is  a  strictly  lower  triangular 
matrix,  and  F  is  a  striotly  upper  triangular  matrix.  With  this 
splitting  of  the  matrix  A,  equation  (2)  can  be  written 

(D  -  E)  x  *  F  x  +  k  (5) 
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Multiplying  both  qldes  of  the  equation  by  a  quantity  CO ,  called  the 
relaxation  factor «  gives 

(COD -COE)  x  =  COF  x  +  CO  k  (6) 

Adding  the  natrlx  product  D  x  to  both  sides  and  rearranging  terns  gives 

<D -COE)  x  ■  £(1  -  (O)  D  +  COfJ  x  +  CO  k  (7) 

This  equation  suggests  the  following  iterative  method: 

(D  -  COE)  x(m41)  =  [(1  -  CO)  £  +  60  f]  x^  +  CO  k,  m  *  0  (8) 

where  the  components  of  x^  are  the  initial  estimates  of  the  unique 
solution  of  x  of  equation  (2).  As  D  -  COE  is  nonsingular  for  any 

■C  4 

choice  of  CO  i  and  with  the  definitions  L  =  D  E  and  |  *D  F,  equation 
(8)  can  be  written 

x(m+l)  -  (I  -COL)"1  [(1  -  CO  )I  +  CO  uj  xU) 

+  CO  ( I  -COL)’1  S"1  k  (9) 

This  Itoratlvo  method!  is  called  the  point  successive  overrelex&tion 

iterative  method  and  will  be  denoted  the  S.O.R.  iterative  method.  The 

matrix 

j*  K  (I  -CO  L)”1  [(1  -  CO  )  I  +  (0  W  J  (10) 

is  called  the  iteration  matrix  of  the  S.O.R.  iterative  method. 

Spectral  Radius.  The  magnitude  of  the  largest  of  the  eigenvalues 
of  the  iteration  matrix  is  termed  the  spectral  radius.  The 
eigenvector  associated  with  this  eigenvalue  will  be  called  the  dominant 
eigenvector  of  the  iteration  matrix.  Varga  (Ref  6:13)  has  shown  that 


% 
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the  S.O.R.  iterative  method  is  convergent  If  and  only  if  the  spectral 
radius 


(11) 


This  condition  la  assured  provided  (a)  that  the  matrix  (D  -  CO  E)  is 

nonsingular  and  its  inverse  is  non-negative,  and  (b)  that  the  matrix 

["(1  -  CO  )D  +  CO  f!  is  non-negative  (Ref  6:89)> 

— (m) 

Error  Vectors.  The  error  vectors  €  associated  tilth  the  S.O.R. 
Iterative  method  are  defined  by 


-(m)  _(m)  _ 

€  =  x  -  x. 


m  &  0 


(12) 


where  x  is  the  unique  vector  solution  of  equation  (2).  By  equation 
(9)  the  error  vectors  can  be  expressed  as 


--(m)  »  — (m-1)  )  a  -(0) 

€  -  Lwe  - - 


(13) 


Assuming  that  the  eigenvectors  of  form  a  complete  set,  the  initial 


=(0) 


error  vector  €  can  be  expressed  in  terms  of  them  as 

_(G) 

€  “  »2«2  +  •  •  •  +  ftn®n 


(14) 


where  the  eigenvectors  ej ,  eg;  •  •  •«  Sp  are  assumed  to  be  normalised. 
The  4i ,  A£r  .  •  .,  an  represent  the  magnitudes  of  the  components  of  the 
initial  error  vector.  Multiplying  aquation  (14)  by  yields 

€U>  -  6<0)  -  *1  +  *2  Lw  ®2  +  •  •  •  *n  «5> 


which  by  the  eigenvalue  equation  e^  =  X  reduces  to 


,(1> 


X  1*1  +  *2  X  2*2  +  •  •  •  +  *n  ^  n®n 


(16) 
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where  X  ^ ,  X  •  •  •  »  X  n  are  the  eigenvalues  of  Continuing 

in  this  fashion  the  error  vector  at  eaoh  iteration  is  given  by 


■eiX2e1+a2X|^2  +  .  .  .+*nX2en 
•  • 

•  • 

•  e 

?<m).MX!['Sl+.2X£i2  +  .  .  .+«n^'*n  (17) 


Assuming  that  the  eigenvalues  are  ordered  by  magnitude  such  that 

|x,|  >  |XJ> •  ••» |x.| 

(18) 

then 

|x:i  >  \k\>  -- 

>  lx:  i 

(19) 

For  sufficiently  large  m,  J  » |  j 

and 

rr(m)  \  at  - 

€  *  «1  A  i  oj. 

(Ref  ifll5) 

(20) 

definition  the  spectral,  radius  p  (  )  is  etjual  to  the  magnitude 

of  the  largest  eigenvalue  X  j  of  L-qj  ,  thus  it  can  be  seen  that  the 
error  vector  will  vanish  as  m  — •»  oo  if  and  only  if  p  (  )  <  i. 

For  sufficiently  large  m  the  magnitude  of  th©  error  vector  at  each 
iteration  is  primarily  a  function  of  the  spectral  radius  and  the 
quantity  a^  defined  by  equation  (14).  Both  of  these  quantities  are 
determined  by  the  choioe  of  the  Iteration  matrix  »  which  in  turn 

Is  determined  by  the  choice  of  the  relaxation  faotor  CJ  and  the  order 
in  which  the  equations  (i)  are  solved*  A  particular  ordering  has  been 
implied  In  this  development  by  the  splitting  of  matrix  A  in  equation 
(4).  This  ordering  is  termed  the  natural  or  standard  ordering. 
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Rate  of  Convergence.  There  is  some  disagreement  in  the  literature 
concerning  the  definition  of  rate  of  convergence.  Young,  in  1950 »  de¬ 
fined  the  rate  of  convergence  as  ~  In  p  (  )  (Ref  7:9*+)»  Varga, 

however,  has  chosen  to  call  this  the  asymptotic  rate  of  convergence 
since  it  is  realized  only  as  in - »  &O  °  Varga  then  defines  the  aver¬ 

age  rate  of  convergence  as 

R(Am)  =  -ln  £  (  1 1  A™  j  j  }1/m  ]  (21) 

for  m  iterations,  where  |  |  A 1 1  denotes  the  spectral  norm  of  A  (Ref  6:62). 

If  two  iteration  matrices,  A  and  B,  have  different  average  rates  of  con¬ 
vergence  for  m  iterations  such  that  R(Am)  <R(rf1)  then  matrix  A  is  said 
to  be  iteratively  faster  for  m  iterations  (Ref  6:62).  Varga's  definitions 
will  be  used  throughout  this  paper.  The  average  rate  of  convergence 
and  the  average  error  reduction  rate  are  equivalent. 

Consistent  Ordering.  Varga  (Ref  6:125)  has  shown  that  any  ordering 
of  equations  (l)  which  is  "consistent"  with  the  natural  ordering  gives 
rise  to  an  iteration  matrix  whose  spectral  radius  is  equal  to  that 

for  h'(jj  £ivsn  in  equation  (10)  for  the  3 are  CiJ  ,  Furthermore,  Varga 
has  shown  that  this  spectral  radius  is  smaller  than  the  spectral  radiuc 
for  any  other  ordering  which  is  not  consistent  with  the  natural  ordering. 
Hence,  for  sufficiently  large  .#  the  rate  of  error  vector  reduction  will 
bo  greatest  for  consistent  orderings. 

Let  the  n  x  n  point  Jacobi  matrix  B  be  defined  by 

B  -  -D'1  A  +  I  (22) 

where  matrix  D  is  defined  in  equation  (4),  and  matrix  A  is  defined  in 
equation  (2).  Then  the  matrix  A  is  consistently  ordered  if  all  the 
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eigenvalues  of  the  matrix 

B(  TJ  )  s  TJ  L  +  TJ  "^p  “  ^  U  (23) 


derived  from  the  matrix  B  =  L  +  U,  where  L  and  U  are  respectively 
strictly  upper  and  strictly  lower  triangular  matrices,  are  independent 
cf  TJ  ,  for  TJ  4  0,  provided  that  A  is  a  p-eyclic  matrix  (Ref  6:101). 
The  matrix  A  is  defined  as  p-cyclie  if  the  matrix  B  of  (22)  is  weakly 
cyclic  of  index  p  (-  2)  (Ref  6:99).  The  matrix  B  is  weakly  cyclic  of 
index  p  (Ref  6:39)  if  there  exists  an  n  x  n  permutation  matrix  P  such 
that  PBPT  is  of  the  form 


The  eigenvalues  of  the  point  successive  overrelaxation  iteration 
matrix  L'qj  derived  from  a  consistently  ordered  p-eyclic  matrix  A 
are  related  to  the  eigenvalues  of  the  associated  Jacobi  matrix  B  by 


<  X  +  CO  -  1>P  =  X^1  COP  jU.P  (25) 

where  \  is  a  nonsero  eigenvalue  of  •  Cl)  is  the  relaxation  factor, 

and  fd.  is  an  eigenvalue  of  B  (Ref  6:106).  That  such  a  relationship 
exists  is  itself  Interesting,  but  its  importance  lies  in  the  fact  that 


knowledge  of  the  eigenvalues  of  the  matrix  B  allows  determination  of 
the  eigenvalues  of .  L'gj  f°r  *  consistent  ordering. 

Directed  Graph.  Fortunately,  it  is  not  necessary  to  apply  the 
definition  directly  to  determine  whether  an  ordering  is  consistent. 
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It  can  be  shown  (Ref  6:121)  that  an 
ordering  Is  consistent  if  it  leads 
to  a  matrix  B  whose  directed  graph 
of  type  2  has  an  equal  number  of 
major  and  minor  paths  on  every 
closed  path*  A  directed  graph  of 
type  2  for  the  matrix  B  =  (b^, j) 
is  constructed  so  that  if  b^j  ^  0 
then  a  path  from  node  P^  to  the 
node  Pj  is  drawn  and  denoted  by  a 
double- arrow  (see  Fig.  1)  only  if  j  >  i;  otherwise,  a  single-arrowed 
path  is  drawn.  The  former  paths  are  called  major  paths;  the  other 
paths  are  called  minor  paths.  It  is  also  possible  to  use  this  graph 
to  determine  when  the  matrix  A  is  p-cyolic;  or  equivalently,  when  the 
point  Jacobi  matrix  B  is  weakly  cyclic  of  index  p  (Ref  6:100).  If  the 
greatest  common  divisor  of  the  lengths  of  closed  paths  of  the  graph  is 
p  then  tho  matrix  B  is  weakly  cyclic  of  index  p,  where  the  length  of  a 
path  is  defined  as  the  number  of  nodes  reached  in  traversing  the  path. 

Optimum  Relaxation  Factor.  As  stated  previously  the  spectral 
radius  of  the  iteration  matrix  depends  on  the  choice  of  the  relax¬ 

ation  factor  CO  ,  as  well  as  the  ordering  of  equations  (1).  Ostrowski 
has  shown  that  the  S.O.R.  iterative  method  is  convergent  for  all  CO 
such  that  0  <  CO  <  2  (Ref  6:77).  The  optimum  CO  ,  in  the  sense  that 
p  (  Lqj  )  is  minimised,  can  be  related  to  the  largest  of  the  eigen¬ 
values  of  the  matrix  B  of  equation  (22),  for  the  special  case  p  =  2,  as 


where  CO  ^  denotes  the  optimum  CO  (Ref  6:110). 


Fig.  1  Type  2 
Directed  Graph 
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Thus,  for  2-cyolie  matrices,  the  problem  of  finding  the  ordering 
and  relaxation  factor  which  yields  the  greatest  rate  of  error  vector 
reduction  in  the  limit  as  m  — *  oo  has  been  reduoed  theoretically  to 
choosing  any  consistent  ordering  and  finding  the  largest  eigenvalue 
of  the  associated  point  Jacobi  matrix  B.  The  theory  still  allows  for 
the  possibility,  however,  that  the  average  rata  of  convergence  may  be 
different  for  different  consistent  orderings,  and  henoe,  that  the 
approach  to  asymptotic  convergence  may  differ  between  consistent  order¬ 
ings.  Also,  the  projection  of  the  initial  error  vector  onto  the 
dominant  eigenvector  of  the  iteration  matrix,  expressed  as  the  quantity 
a^  in  equation  (20),  may  be  different  for  different  orderings  even  if 
all  are  consistent.  The  extent  to  which  this  may  be  of  practical 
importance  has  been  examined  In  this  study  by  applying  the  point  suc¬ 
cessive  overrelaxation  iterative  method  to  the  finite  difference 
approximations  to  elliptic  and  parabolic  partial  differential  equations. 

Difference  Equations 

To  derive  these  difference  equations  consider  the  general  second 
order,  linear,  partial  differential  equation 

A  02  u  +  B  0  2  u  +  C  02  u  +  D  0  u  +  E  0  u  +  Fu  =  G 

0X2  0X0  y  0  y^  0  x  0  y  (27) 

where  A,  B,  ...  ,  G  are  constants  or  function  of  x  and  y  only.  This 
equation  is  classified  as  elliptic,  parabolic,  or  hyperbolic  in  a 
domain  of  the  xy  plane  as  the  values  of  the  function  B2  -  4  A  C  are 
negative,  zero,  or  positive,  respectfully,  throughout  the  domain. 
Equations  of  the  hyperbolic  type  were  not  considered  in  this  study. 


% 
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Laplace’s  Equation*  Laplace’s  equation  is  an  exunple  of  the 
elliptic  type,  and  this  equation  was  chosen  for  study  bocause  the 
eigenvalues  of  the  Iteration  matrix  of  the  associated  difference 
equations  can  be  found  theoretically.  This  allows  precise  determina¬ 
tion  of  the  optimum  overrelaxation  factor  prior  to  aotual  solution  of 
the  difference  equations. 

To  derive  a  finite  difference  approximation  to  Laplace's  equation 

02  u  (x,y)  +  0^  u  (x,y)  =  o 

0  x2  0  y2  (28) 

subdivide  the  x-y  plane  into  sets  of  equal  rectangles  of  sides 
S  x  »  h,  S  y  =  k,  as  shown  in  Fig.  2,  and  let  the  coordinates  (x,y) 
of  a  representative  mesh  point  p  be 

x  »  ih,  y  **  jk  (29) 

where  i  and  J  are  Integers  (Ref  5:?)* 


y 

Jk 

1 

i,  j  ♦  1 

,p(ih.;IO  . 

1  .  .  J 

<.J 

T 

k 

± 

0 

1* - h - * 

i  h  x 

Fig.  2  Sample  Mesh 
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Denote  the  value  of  u  at  p  by 


^  =  u(  jk)=  u^j  (30  ] 

Expand  the  function  u(x,y)  in  the  x  direction  about  the  point  x^  in  a 
Taylor's  expansion  as 


<HXp  +  h,yp)  = 


u(xp,yp)  +  h  S.  u  +|h2  c)2  u 

0  x  *p.yp  0  x2  Xp.yp 

+  1  h3  03  ^ 

—  —  1  ~  -  "  T  •  •  • 

6  0  ry  Xp,yp 


Addition  of  the  expansions  for  u(xp  +  h9yp)  and  u(xp  -  h,yp)  gives 


u(Xp  +  h,yp)  +  u(xp  -  h,yp)  =  2u(xp,yp)  +  h2  02  « 

a*2  vyP 

+  0(h4) 


where  0(h4)  denotes  terras  containing  fourth  and  higher  powers  of  h. 
Assuming  that  these  terras  are  negligible  compared  with  the  lower 
powers  of  h,  it  follows  that 

d  2  u  -  u(Xp  +  h9yp)  -  2u(xp,yp)  +  u(xp  -  htyp) 

0X2  Xp,yp  h2  (33) 

In  terras  of  the  notation  introduced  in  equation  (29)  this  can  be 
expressed  as 

0 2  u  _u[(i  +  l)h,jk]  -  2u  ih, jkl  +  u[(i  -  Dh.jkl 

0x2  p  h2  (34) 


0*2  i.J 


=  ^*1.-1  ~  2ui.1  *  Uj-1, 
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The  derivation  for  0  2  u/  0  y^  is  similar  and  gives 
d2  u  =  «i, j+1  -  2ui. J  +  “i.J-l 

0  y2  A» J  ic2  06) 

Both  equations  (35)  and  (36)  are  approximations  to  the  derivatives , 
and  are  called  finite  difference  approximations.  The  magnitude  of 
the  error  introduced  by  these  equations  is  reduced  as  the  mesh  size 
is  deoreased. 

Substituting  the  expressions  for  the  derivatives,  equations  (35) 
and  (36),  gives  the  following  approximation  to  equation  (28),  Laplace's 
equation. 


“i+l.j  “  +  “i-l.J  +  **1,3+1  “  +  ui»j-l  *  0 


Assuming  for  simplicity  of  analysis  that  the  mesh  spaclngs  are 
taken  equal  in  the  x  and  y  directions,  that  is,  that  h  =  k,  equation 
(37)  reduces  to 

4ui.  j  ~  <ui+l,j  +  «i-l,j  +  "i,  j*l  +  "i.j-l)  B  0  (38) 

for  each  mesh  point.  This  expression  is  called  a  five-point  formula 
since  each  mesh  point  (xj.,yj)  is  coupled  to  at  most  four  other  adja¬ 
cent  mesh  points.  This  system  of  linear  equations  can  be  written 

A  u  =  k  (39) 

where  the  column  vector  k  oontains  the  values  of  the  function  u 
specified  on  the  boundaries  of  a  rectangular  region  in  the  xy  plane. 
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The  derivation  leading  to  equation  (39)  shows  that  A  is  a  real 
symmetric  matrix  with  positive  diagonal  entries  and  non-positive  off- 
diagonal  entries.  Moreover,  it  oan  be  shown  that  A  is  irreduoibly 
diagonally  dominant,  so  that  A  is  positive  definite  (Kef  6:18?). 

This  implies  that  A”*  >  0,  where  £  is  the  null  matrix  (Ref  6:85). 

The  associated  point  Jacobi  matrix  B,  defined  by  equation  (22), 
oan  be  shown  to  be  a  non-negative.  Irreducible  cyclic  matrix  of  index 
2  with  real  eigenvalues,  and  with  p  (B)  <1  (Ref  6:188). 

Since  all  of  the  diagonal  entries  of  the  matrix  A  are  4,  the 
point  Jacobi  matrix  B  associated  with  A  can  be  expressed  simply  as 

B  =  I  -  iA  (40) 

All  tke  eigenvalues  and  eigenvectors  of  the  matrix  B  can  be  obtained 
explicitly  for  the  finite  difference  approximation  to  Laplace's 
equation  in  a  rectangle.  Since  B  is  non-negative  and  irreducible, 
it  has  a  unique  positive  eigenvector  with  eigenvalue  p  (B) ,  which 
for  a  unit  square  is  given  by 

p(B)  =  cos  (h7T)  (Ref  6:203)  (41) 

For  the  point  successive  overrelaxation  iterative  method  the 
optimum  CO  can  thus  be  written,  by  equation  (26),  as 

,  *  2  2 

COjj  =  - -  - 

l+\/l-p2(B)  1  +  sin  (hTT)  (42) 

for  all  consistent  orderings. 

Heat  Equation.  The  derivation  of  the  five  point  finite  differ¬ 
ence  approximation  to  the  diffusion  equation  of  transient  heat  transfer 
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is  similar  to  that  for  Laplace's  equation.  The  mesh  s pacings  and 
notation  defined  in  equation  (30)  and  illustrated  in  Fig.  2  are  again 
adopted. 

In  two  dimensions  the  heat  equation  for  an  isotropic,  homo¬ 
geneous  material  can  be  written  as 


0  U  =  <Y  0  2  0  +  02  U 
0  t  0  07 


(43) 


where 

U(x,y,t)  is  the  temperature, 

CL  is  the  thermal  diffusivity,  and 

depends  an  the  rate  of  heat  generation,  the  density  of 
the  material,  and  the  heat  oapaoity. 

It  can  be  seen  that  the  steady  state  heat  equation  with  no  heat  gener¬ 
ation  reduces  to  Laplace's  equation,  while  the  transient  heat  equation 
is  parabolic. 

With  the  assumption  ef  niform  mesh  spacing,  equal  to  h,  the 
second  partial  s  of  0  with  respect  to  x  and  y  at  time  t  +  At  can  be 
replaced  with  finite  differences  similar  to  equations  (35)  end  (36) 
while  the  first  partial  of  0  with  respect  to  t  can  be  found  from  the 
Taylor's  expansion 


UCxp.yp.t)  ■  U(xp,yp,t  ♦  A  t)  -  At  0P 

0t 


Xp»yp»t  ^  A  t 


+  i(  At)2  0  2  u 


xp.yp.t  +  A  t 


(44) 


Assuming  that  higher  powers  of  A  t  can  be  neglected,  this  becomes 
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a t 


Xp.yp.t  +  A  t 


(45) 


Using  these  finite  difference  equations  at  tine  t  +  At,  equation  (43) 
can  be  expressed  as 


Introducing  the  notation  r  3  CL  At/h^  this  equation  becomes 

(1  +  4r)  0(i,j,t+At)  *  r  £u(i«*l,j,t+ At)  +  D(i-l,J,t+At) 

+  D(i,j44ft+At)  +  U(i,j-l,t+At)] 

+  U(ifj,t)  +  jS  At  (47) 


or 

U(i,j,t+ At)  «_r_  P U(l+1  f jft+  At)  +  0(W,J,t*At) 

l"*4r 

+  u(i,j-n,t+At)  +  o(i,j-i.t+At)] 

+  _L_  fua.^t) +/?  At] 

144r  L  ^  J  (48) 


This  system  of  linear  equations  can  be  written 

A  U  =  K  (49) 

where  K  includes  the  given  boundary  and  Initial  conditions  as  well  as 
the  heat  generation  rate.  Again,  as  for  Laplace's  equation,  A  is 
irreducibly  diagonally  dominant,  positive  definite,  and  thus  A“*  >  0. 
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The  associated  point  Jaoobi  matrix  B  is  a  non-negative,  irreducible 
cyclic  matrix  of  index  2  with  real  eigenvalues,  and  with  p  (B)  <1. 

The  eigenvalues  of  the  matrix  B,  however,  cannot  be  determined 
explicitly  for  this  equation*  Thus-  the  spectral  radius  O  (B)  must 
be  estimated,  or  the  optimum  GO  determined  experimentally,  for  a 
given  problem. 

Ordering 

The  standard  method  of  ordering  of  the  difference  equations, 
such  as  equation  (38),  for  solution  by  an  iterative  method  can  be 
illustrated  for  a  small  sample  mesh,  as  shown  in  Fig.  3*  The  node 
(i  *  1,  J  *  1)  is  assigned  the  first  position  in  the  ordering,  the 
node  (1,2)  the  second,  and 
(1,3)  the  third.  The 
"scan”  is  then  continued 
in  a  similar  manner  for 
the  row  i  =  2  giving  node 
(2,1)  the  fourth  position 
in  the  ordering,  node  (2,2) 
the  fifth  and  so  on 
throughout  the  mesh.  If 
each  node  within  the  mesh 
is  assigned  a  number  as 
shown  in  Fig.  3,  the  standard  ordering  may  be  expressed  as  ^1,2,3, 
^•5*6,7, 8,9 J  ,  where  the  first  position  within  the  brackets  design¬ 
ates  which  point  is  to  be  solved  first,  the  second  position  designates 
which  is  to  be  solved  second  and  so  on. 


*-h— ► 
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h 

i 

7 

8 

9 

A 

5 

6 

1 

2 

3 

O  12  3  4 
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Fig.  3  Mesh  Numbering 
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With  this  notation  *  different  ordering (  called  the  even- odd 
scan*  can  be  expressed  ^2,4,6,8,1,3»5»7»9 J  •  Again,  the  numbers 
within  the  brackets  denote  the  points  as  labled  in  Fig.  3,  while  the 
positions  within  the  brackets  denote  the  order  in  which  the  points 
are  taken  for  solution,  with  the  first  position  first,  the  second, 
second,  and  so  on.  When  the  point  occupying  the  last  position  is 
solved  for,  the  iteration  is  complete.  The  ordering  is  also  called 
the  "scan". 

The  definition  of  a  consistent  ordering  can  new  be  illustrated 
for  the  two  scans  defined  above.  Take,  for  aocanple,  the  finite  diff¬ 
erence  approximation  to  Laplace's  equation  in  a  square.  The  mesh 
and  point  numberings  shown  in  Fig.  3  will  be  used.  Equation  (38)  can 
be  written  in  matrix  notation  as 

A  u  =  k  (50) 


where  the  column  vector  k  contains  the  values  of  the  funotion  u 
specified  on  the  boundaries  of  the  mesh.  The  matrix  A  is  given  by 


A  = 


l  -£  o  -i  o 

l  — b  o  -£ 
o  -i  l  o  o 

-i  o  o  l  -i 

0  — ^  0  — ^  1 

o  o  -i  o  -i 

o  o  o  -i  o 

o  o  o  o  -i 

ooooo 


0  0  0  0 

0  0  0  0 

-i  o  o  o 

o  -i  o  o 

-i  o  -i  o 

i  o  o  -i 

o  l  -i  o 

0  — ^  1 

-i  o  -i  l 


(5D 
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The  associated  Jacobi  matrix  B  is  defined  by 

B  =  1  -  D*1  A 


(52) 


and  is  given  by 


(53) 


To  determine  whether  an  ordering  is  consistent,  the  matrix  B  is 

transformed  by  the  permutation  matrix  P  corresponding  to  the  ordering, 

-1  T 

by  a  similarity  transformation  of  the  form  P BP  =  PBP  (Ref  3:2h4). 

The  permutation  matrix  P  corresponding  to  the  standard  ordering 

£l,2, 3,  *r,5, 6, 7, 8, 9]  is  simply  P  =  I  and  hence  the  transformation 

yields  the  matrix  B.  Now,  the  ordering  is  consistent  if  the  directed 

T 

graph  of  type  2  of  the  matrix  PBP  has  an  equal  number  of  major  and 
minor  paths.  The  directed  graph  of  type  2  for  matrix  B  of  equation 
(53)  is  shown  in  Fig,  4, 

Many  useful  properties  can  be  deduced  from  this  graph.  The  graph 
is  "strongly  connected",  that  is,  there  is  a  path  (with  arrows)  from 
each  point  to  every  other  point.  Hence  the  matrix  B  (and  therefore  A) 
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Fig.  4  Directed  Graph  of  Type  2 
for  the  Standard  Scan 


Is  irreducible  (Ref  6:20).  The  greatest  common  divisor  of  closed  paths 

is  2,  hence  the  matrix  B  is  weakly  cyclic  of  index  2  (Ref  6:101).  By 

definition,  then,  the  matrix  A  is  2-cyelio  (Ref  6:99)»  The  number  of 

major  (2- arrowed)  and  minor  (1 -arrowed)  paths  on  every  closed  path  is 

equal;  therefore,  the  matrix  B  is  consistently  ordered.  Since  the  same 

T 

conclusions  hold  for  PBP  ,  the  standard  scan  is  consistent. 

The  permutation  matrix  P  for  the  even-odd  scan  £ 2, 4, 6, 8,1 ,3»5»7»9j 
is  given  by 


010000000 

000100000 

000001000 

000000010 

P=  100000000  (54) 

OOICOOOOO 
000010000 
000000100 
_0  0  0  0  0  0  0  0  1_ 

T 

The  directed  graph  of  type  2  for  PBP  .  where  3  is  given  by  equation 
(53),  is  shown  in  Fig.  5«  The  number  of  major  and  minor  paths  on  each 
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Fig.  5  Directed  Graph  of  Type  2 
for  the  Even- Odd  Sean 


closed  path  Is  equal,  and  therefore  the  even- odd  scan  Is  consistent. 
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III..  Experimental  Investigation  and  Results 

In  this  investigation  the  standard  scan  and  the  even-odd  scan 
were  programmed  in  Fortran  IV  for  execution  on  an  IBM  ?09^  computer. 

The  scans  were  applied  to  the  five  point  finite  difference  approxima¬ 
tions  to  Laplace's  equation  and  the  heat  equation  in  a  square  with  a 
uniform  mesh  spacing  in  both  directions.  The  method  of  point 
successive  overrelaxation  was  used. 

These  two  scans  were  compared  to  each  other  and  to  a  method  devel¬ 
oped  and  tested  by  Cudahy  in  1965  (Ref  2).  Cudahy  applied  his  method, 
which  was  designed  to  accelerate  the  convergence  of  the  standard  point 
successive  overrelaxation  iterative  method,  to  the  five  point  finite 
difference  approximation  to  the  diffusion  equation  of  transient  heat 
transfer,  using  point  successive  overrelaxation.  He  compared  his  method 
to  the  standard  scan  and  found  that  his  method  was  iteratively  faster, 
the  per  cent  difference  depending  on  various  factors,  including  the 
time  step  and  heat  generation  rate.  lie  did  not  program  the  even-odd 
scan. 

One  of  the  first  objectives  of  the  present  study  was  to  explain, 
from  an  examination  of  the  theory  of  point  successive  overrelaxation, 
why  Cudahy's  method  was  faster  than  the  standard  scan.  This  proved 
to  be  somewhat  difficult,  however,  since  the  theory  that  has  been 
developed  for  point  successive  overrelaxation  deals  with  cyclic, 
stationary  iterations,  while  Cudahy's  method  strictly  speaking  is 
neither  cyclic  nor  stationary. 

An  iterative  proced\ire  is  said  to  be  cyclic  when  the  ordering  is 
repeated,  without  change,  on  each  iteration.  The  method  developed  by 
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Cudahy  employs  four  different  orderings  (each  of  which,  considered  by 
itself,  is  consistent).  These  four  scans  are  applied  repetitively,  in 
the  same  order,  until  the  process  converge  -  to  a  predetermined  limit. 
The  four  orderings  can  be  illustrated  in  terms  of  the  mesh  points 
labeled  in  Fig.  3.  The  fir/t  ordering  is  £l,2,3»7»8,9»^»5*6j  ;  the 
second  ordering  is  |^3»2,1 ,9(8, 7,6, ;  the  third  ordering  is 
[l»^.7»3t6,9.2,5.8]  ;  the  fourth  ordering  is  [^7, 4, 1,9, 6, 3*8,5, 2  J  . 
The  center  row  or  column  is  always  solved  last. 

An  iterative  process  is  termed  stationary  if  the  relaxation 
factor  applied  is  independent  of  the  iteration  number.  In  his  method, 
Cudahy  applied  two  different  relaxation  faotors  on  each  iteration. 
Although  the  two  relaxation  factors  are  the  same  on  each  subsequent 
iteration,  the  use  of  four  different  orderings  makes  the  relaxation 
factor  applied  to  certain  points  of  the  mesh  a  function  of  the 
iteration  number,  and  hence  the  method  is  nan-stationary. 

One  of  the  difficulties  that  this  presents  oan  be  seen  by 
examining  the  assumptions  used  in  the  error  vector  analysis  leading  to 
equation  (20 )  of  chapter  II.  It  was  assumed  there  that  the  same  itera¬ 
tion  matrix  L of  equation  (10)  was  applied  on  each  iteration,  that 
is,  that  the  method  was  cyclic  and  stationary.  For  Cudahy's  method 
there  are  four  different  iteration  matrices,  which  do  not  have  the  same 
set  of  eigenvectors,  and  thus  equation  (20)  does  not  apply. 

Another  difficulty  introduced  by  the  use  of  four  iteration 
matrices  is  encountered  in  predicting  the  asymptotic  rate  of  conver¬ 
gence.  Tho  rate  is  affected  by  a  coupling  between  the  iteration 
matrices  which  cannot  be  exactly  known  without  knowledge  of  the  full 
range  of  eigenvalues  and  eigenvectors  of  each  laairix.  Determination 
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of  the  eigenvalues  nd  eigenvectors  of  the  iteration  matrices,  which 
is  a  formidable  p.  o-*-sm  in  itself,  is  complicated  by  tlv  fact  that 
eocplicit  representation  of  i  e  iteration  matrices  requires,  for  a  51 
by  51  mesh,  the  direct  inversion  of  a  2401  by  2401  matrix  for  each 
iteration  matrix.  This  type  of  problem  is  avoided  in  the  theory  of 
point  successive  overrelaxation  by  equation  (25).  which  specifies  the 
eigenvalues  of  the  S.O.R.  iteration  matrix  in  terns  of  the  eigenvalues 
of  the  associated  Jacobi  matrix  B,  which  can  be  explicitly  represented 
from  equation  (22). 

Since  the  theory  provides  no  immediate  answer  for  these  problems, 
the  computer  was  used  to  compare  Cudahy’s  method  experimentally  with 
the  standard  and  even-odd  scans. 

Laplace’s  Equation 

The  finite  difference  approximation  to  Laplace's  equation,  for  a 
uniform  mesh  spacing,  given  by  equation  (38),  can  be  solved  by  the 
method  of  point  successive  overrelaxation  derived  in  chapter  II.  The 
point  equation  corresponding  to  1  1  matrix  equation  (9)  is  given  by 

u(i, j)  =  (1  -  CO  )u(i,j)  +  CO  [~u(l  +  l,j)  +  u(i  -  l,j) 

4  L 

+  u(i, j  +  1)  +  u(i, j  -  1)J  (55) 

where  the  (m  +  l)'s  and  (m)'s  have  been  suppressed  to  allow  this 
equation  to  represent  any  ordering.  The  values  of  the  fun  .on  u 
within  the  brackets  are  to  be  taken  from  the  present  or  preceding 
iteration,  depending  on  the  ordering  of  points  for  solution. 

Equation  (55)  is  applied  at  each  point  within  the  mesh  to  form  a 
complete  iteration.  The  magnitude  of  the  ohange  in  the  value  of  u  at 
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each  point  Is  **  ted  as  the  Iteration  progresses*  and  the  largest 
change  Is  retained.  At  the  end  of  each  iteration  this  value  can  be 
compared  to  a  predetermined  number*  called  the  Iteration  limit.  If 
the  largest  of  the  changes  In  magnitude  of  the  function  u*  denoted  y  * 
Is  less  than  the  iteration  limit  the  process  may  be  said  to  have  con¬ 
verged  to  that  iteration  limit.  This  is  an  easily  applied  convergence 
criteria. 

In  this  study  the  difference  magnitude  y  was  not  tested  at  the 
end  of  each  iteration  to  determine  convergence*  but  rather  a  specified 
number  of  iterations  were  performed  with  the  value  of  y  printed  out 
for  eaoh  iteration.  If  an  iteration  limit  is  later  specified  then  the 
first  iteration  where  y  is  less  than  the  specified  limit  will  be  the 
Iteration  at  which  the  process  can  be  said  to  have  converged. 

The  standard  and  even- odd  scans  where  compared  to  Cudahy’s  method 
for  a  31  by  31  mesh  where  the  value  u  =  1000  was  specified  on  the 
boundaries,  and  the  initial  estimate  for  all  points  within  the  mesh 
was  u  *  0.  The  results  of  tha<  omparison  are  shown  in  Table  1  for 
various  iteration  limits.  It  was  stated  previously  that  Cudahy's 
method  employs  two  relaxation  factors.  The  second  relaxation  factor 
Is  applied  only  to  the  points  In  the  last  row  or  column  of  the  iteration 
and  is  given  in  terms  of  the  O  listed  in  Table  1  by 

60 ’  *  frX1  +  (Ref  2:25) 

COr  +  1  +  4r  (56) 

where  r  =  Q,  ^  which  gives  for  Laplace's  equation  (At - *  00) 

h2 

Cl)'  "  4Cl> 

4  +CJ  (57) 
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Table  i 


Number  of  Iterations  Required  for  Various  Iteration 
Limits  for  a  31  x  31  Mesh  Solving  Laplace's  Equation 


Iteration 

Limit 

Iterations 

Standard 

Soar. 

Even-Odd 

Scan 

Cudahy's  Method 

60  =  1.81 

00 

• 

II 

3 

60*1.80 

.  J 

60*1.85 

60*1.90 

60*1.95 

102 

18 

16 

5 

5 

9 

24 

101 

36 

28 

31 

27 

31 

67 

H. 

o 

o 

61 

40 

65 

55 

64 

94 

10”1 

62 

y 

99 

82 

75 

137 

10”2 

71 

64 

133 

110 

94 

191 

10“3 

85 

76 

167 

138 

128 

232 

10"4 

9  6 

88 

202 

166 

144 

288 
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It  can  be  seen  from  Table  1  that  Cudahy’s  method  is  initially  itera¬ 
tively  faster  than  either  of  the  other  two  scans.  As  the  iteration 
limit  becomes  finer,  however,  Cudahy's  method  becomes  iteratively 
slower.  Thus  it  is  likely  that  the  asymptotic  rate  of  convergence  of 
Cudahy's  method  is  less  than  that  for  a  consistent,  cyclic  stationary 
iteration.  This  assertion  cannot  be  proven  experimentally,  of  oourse, 
since  the  asymptotic  rate  of  convergence  is  realized  only  as  the 
number  of  iterations  approaches  infinity. 

Since  the  standard  scan  and  the  even-odd  scan  are  consistent,  the 
optimum  overrelaxation  faotor  can  be  determined  precisely  from  equation 
(42).  The  optimum  CO  was  used  to  obtain  the  data  shown  in  Table  1  for 
these  scans.  It  is  interesting  to  note  that  the  even-odd  scan  is  iter¬ 
atively  faster  than  the  stanuard  scan  for  all  iteration  limits  shown. 
Since  both  scans  have  the  same  asymptotic  rate  of  convergence,  the  even- 
odd  scan  can  be  expected  to  be  iteratively  faster  than  the  standard  scan 
for  all  iteration  limits,  no  matter  how  fine. 

Figs.  6  and  7  show  a  segment  of  the  data  from  which  Table  1  was 
constructed,  for  the  standard  and  even-odd  scans.  The  difference 
magnitude  'Y  ,  defined  as 

y  E  max  J  u<“*>  -  |  ,  1  -  i,  j  -  •  t  (58) 

where  N  is  the  number  of  mesh  spacings,  is  plotted  versus  the  Iteration 
number.  For  a  sufficiently  large  number  of  iterations  the  even-odd  scan 
is  iteratively  faster  than  the  standard  scan  by  an  amount  which  is  repre¬ 
sented  by  a  fixed  number  of  iterations  at  a  convergence  rate  approx¬ 
imately  equal  to  the  asymptotic  rate  of  convergence.  Since  the 
difference  between  the  scans  is  a  fixed  number  of  iterations,  the  per 
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Iteration  Number 


Fig*  6  Difference  magnitude  Y  versus  iteration  number 
_ for  a  31  x  31  mesh  solving  Laplace^  equation 
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Fig.  7  Difference  magnitude  'Y  ▼eraua  iteration  number 
for  a  31  *  31  mesh  solving  Laplace's  equation 
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cent  difference  is  a  function  of  the  number  of  iterations  and  will 
approach  zero  as  the  number  of  iterations  approaches  infinity. 

Fig.  8  shows  a  comparison  between  the  even-odd  scan  at  optimum 
CO  and  Cudahy's  method  at  (O  =  1.85.  It  can  be  seen  that  Cudahy's 
method  rapidly  approaches  its  asymptotic  rate  of  convergence,  which 
is  not  as  great  as  that  for  the  even-odd  scan.  The  somewhat  erratic 
behavior  of  the  difference  magnitude  for  the  first  few  iterations 
illustrates  the  difficulty  in  making  theoretical  predictions  of  aver¬ 
age  rates  of  convergence  for  a  relatively  small  number  of  iterations. 

To  determine  whether  the  mesh  size  has  an  important  effect  on  the 
relative  speed  between  the  standard  and  even- odd  scans,  the  mesh  size 
was  varied  from  a  10  by  10  to  a  70  by  70.  Computer  time  limitations 
precluded  inclusion  of  larger  mesh  sizes.  The  results  of  that  compar¬ 
ison  are  shown  in  Table  2.  Again  the  optimum  CO  was  calculated  from 
equation  (42)  and  used  in  these  runs.  The  iteration  limit  used  repre¬ 
sents  a  reduction  in  the  length  of  the  error  vector  by  a  factor  of 
approximately  10”**,  or  eleven  orders  of  magnitude. 

The  data  in  Table  2  show  that  the  number  of  iterations  difference 
between  the  two  scans  increases  as  the  mesh  spacing  decreases,  and  that 
the  per  cent  difference,  for  the  same  iteration  limit,  remains  approx¬ 
imately  the  same.  This  indicates  that  the  relative  speed  difference 
is  independent  of  mesh  spacing. 

Figs.  9  and  10  show  the  difference  magnitude  y  as  a  function  of 
the  iteration  number  for  a  40  by  40  mesh.  It  can  be  seen  that  the 
same  general  features  are  present  here  as  were  present  for  the  31  by  31 
mesh  in  Figs.  6  and  7.  The  only  significant  change  is  the  difference 
in  asymptotic  rate  of  convergence,  the  slope  of  the  curve  for  large 
iteration  numbers. 
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Table  2 


Number  of  Iterations  Required  at  an  Iteration 
Limit  of  10“°  for  Solution  of  Laplace’s  Equation 
with  Various  Mesh  Sizes  for  a  Unit  Square 


I 


Mesh 

Size 

Iterations 

Standard 

Scan 

Even-Odd 

Scan 

Difference 

Difference 

10  x  10 

48 

4  5 

3 

6.3# 

20  x  20 

98 

93 

5 

5.1# 

30  x  30 

148 

139 

9 

6.1* 

40  x  40 

195 

186 

9 

4.656 

50  x  50 

244 

231 

13 

5.3# 

60  x  60 

293 

2  77 

16 

5.4# 

70  x  70 

343 

323 

20 

5.8# 
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Fig.  10  Difference  magnitude  'Y  versus  iteration  n  rmbor 
for  a  40  x  40  mesh  solving  Laplace's  equation 
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Heat  Equation 

The  point  equation  for  the  finite  difference  approximation  to 
the  heat  equation,  corresponding  to  the  matrix  equation  (9),  is  given 

t>y 

u(i,j,t+  At)  ^  (i-co)u(i,j,t+ At)  +  a)  |  j  [u(i+i,j,t+At) 

\  l"*4r 

+  u(i-i,j,t+ At)  +  u(i,.tu,t+At)  +  u(i,>i,t+At)] 

+  1  [o(i,j,t)  +  ft  At]  \ 

l+t*  '  j  (59) 

where  the  (ra  +  l)'s  and  (m)'s  have  been  suppressed  to  allow  this 
equation  to  represent  any  ordering. 

Solutions  of  equation  (59)  by  point  sue  'Hive  overrelaxation  were 
obtained  using  the  standard  and  even-odd  scans  and  Cudahy's  method. 

The  problem  selected  was  similar  to  that  for  Laplace's  equation.  The 
boundary  values  were  set  at  0  =  1000,  the  thermal  diffusivity  Q,  was 
taken  as  Q,  =  1,  and  there  was  no  heat  generation  (  j3  s  0).  This 
corresponds  to  the  problem  for  which  Cudahy  reported  the  greatest 
improvement  in  rat*  of  convergence  for  his  method  compared  to  the 
standard  scan.  The  initial  values  inside  the  mesh  were  0=0  for  each 
mesh  point.  A  uniformly  spaced  51  by  51  mesh  was  employed. 

Tables  3»  4,  and  5  show  the  number  of  iterations  required  for 
various  iteration  limits  foi  At  =  0.01,  0.1,  and  0.5  respectively. 
These  results  correspond  closely  with  those  given  by  Cudahy  (Ref  2) 
for  his  method  and  the  standard  scan  up  to  an  iteration  limit  of  10"”*. 
The  present  data  have  been  tended  to  much  finer  iteration  limits, 
however,  in  order  “o  indicate  the  asymptotic  behavior  of  the  methods. 

The  even-odd  scan,  which  was  founP  o  be  iteratively  faster  than  the 


36 


GSP/ph/68-5 


Table  3 


Number  of  Iterations  Required  for  Various  Iteration  Limits  for 
a  51  x  51  Mesh  Solving  the  Heat  Equation  with  At  =  0.01 


Iteration 

Limit 

Iterations  j 

Standard 

Scan 

Even-Odd 

Scan 

Cudahy's 

Method 

CO  =  1.25 

CO  =  1.25 

CO  =  1.35 

102 

3 

3 

3 

101 

6 

5 

4 

0 

0 

ri 

9 

7 

6 

1 

O 

12 

8 

8 

CM 

1 

O 

13 

10 

10 

C'* 

1 

0 

18 

11 

12 

10-4 

21 

13 

14 

10”5 

25 

15 

16 

10“6 

28 

16 

18 

y-L 

O 

1 

31 

18 

21 

O 

1 

00 

35 

20 

23 

10"9 

38 

21 

25 

0 

1 

0 

41 

23 

28 

H 

1 

O 

H 

44 

25 

31 

10"12 

48 

26 

33 
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Table  4 


number  of  Iterations  Required  for  Various  Iteration  Limits  for 
a  51  x  51  Mesh  Solving  the  Heat  Equation  with  ^t  =  0.1 


Iterations 

Iteration 

Limit 

Standard 

Scan 

Even-Odd 

Sc  n 

Cudahy’s  Method 

Cl)  =1.64 

C 0  ~  .64 

Cl)  =1.68 

CO  =1.70 

co=i .  72 

102 

9 

7 

4 

4 

4 

101 

18 

12 

7 

7 

7 

0 

0 

tH 

28 

17 

12 

11 

11 

10_1 

38 

22 

18 

17 

16 

10“2 

48 

27 

25 

23 

20 

10"3 

58 

32 

32 

29 

27 

-Jl 

10^ 

68 

37 

39 

35 

32 

1 

O 

▼H 

78 

42 

4  7 

42 

41 

10~6 

89 

4? 

55 

50 

53 

• 

0 

.-t 

100 

53 

64 

60 

66 

00 

l 

0 

T-l 

102 

59 

73 

70 

79 

10'9 

105 

63 

81 

82 

95 

10~10 

109 

68 

90 

96 

108 

10**11 

112 

73 

99 

108 

115 

w 

•rl 

1 

O 

116 

78 

111 

115 

123 
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Table  5 

Number  of  Iterations  Required  for  Various  Iteration  Limits  for 
a  51  x  51  Mesh  Solving  the  Heat  Equation  with  At  =  0.5 


Iterations 

•« 

Iteration 

Limit 

Standard 

Scan 

Even-Odd 

Scan 

Cudahy's 

Method 

00  =  1.79 

00  =  1.79 

00  =  1.85 

102 

16 

12 

5 

101 

34 

22 

11 

10° 

52 

31 

23 

rl 

1 

O 

69 

43 

41 

10“2 

90 

55 

61 

10“3 

101 

67 

90 

10-* 

104 

79 

108 

f 

O 

T-i 

112 

91 

117 

10-6 

123 

103 

131 

10-7 

135 

115 

144 

00 

t 

0 

H 

148 

127 

160 

10-9 

160 

_ 

140 

178 
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standard  scan  for  Laplace's  equation,  is  seen  here  to  be  iteratively 
faster  for  the  heat  equation  as  veil. 

Again,  as  was  the  case  for  Laplace's  equation,  Cudahy's  method 
is  initially  iteratively  faster  than  either  of  the  two  consistent, 
cyclic,  stationary  scans.  The  introduction  of  the  factor  At  in  the 
heat  equation  has  the  effect  of  expanding  the  number  of  Iterations, 
in  relation  to  a  fixed  convergence  criteria,  for  which  Cudahy's  method 
is  faster  than  the  standard  soan.  Comparing  the  plots  for  the  standard 
scan  in  Figs.  7  and  IQ  with  the  data  presented  in  Table  4,  it  can  be 
seen  that  the  sharp  drop  in  the  difference  magnitude  'y ,  which  for 
Laplace's  equation  occur ed  at  an  iteration  limit  of  approximately  10^, 
occurs  for  a  At  =  0.1  at  an  iteration  limit  of  approximately  10“^. 
Cudahy's  method  is  faster  than  the  standard  soan  so  long  as  comparisons 
are  confined  to  this  region,  and  not  extended  to  finer  iteration  limits. 
The  even-odd  scan,  on  the  other  hand,  is  not  adversely  affected  by  the 
transition  to  a  small  At. 
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IV.  Discussion  and  Conclusions 


Consistent  Scans 

As  shown  by  the  tables  and  figures  in  chapter  III,  the  even-odd 
scan  was  found  to  be  iteratively  faster  than  the  standard  scan  for 
all  problems  tested.  This  effect  can  be  explained  in  terms  of  the 
theory  presented  in  chapter  II  by  a  consideration  of  the  expansion  of 
the  error  vector  at  the  (m)th  iteration  in  terms  of  the  eigenvalues 
and  eigenvectors  of  the  iteration  matrix. 

For  a  sufficiently  large  numoer  of  iterations  the  error  vector 
following  the  (ra)th  iteration  is  given  by  equation  (20)  as 

— (m)  \  m  —  ,  . 

€  =  a^  A  i  ej  (20) 

This  condition  can  be  found  to  be  satisfied  experimentally  when  the 
ratio  of  the  difference  magnitude  'y  for  +he  (m)th  iteration  to  the 
difference  magnitude  'y  for  the  (m  -  1  )t.,  ‘■•oration  is  approximately 
equal  to  the  spectral  radius  and  when  this  value  is  essentially  con¬ 
stant  for  several  preceding  and  several  following  iterations.  For  the 
problems  tested  this  condition  was  found  to  apply  for  all  iterations 
past  the  100th  iteration.  Expressed  another  way,  X  was  found  to 
be  much  less  than  X 

In  examining  the  terns  in  equation  (20)  it  will  be  recalled  that 

the  dominant  eigenvector  e^  was  assumed  to  be  normalized.  The  largest 

eigenvalue  X  ^  is  the  same,  and  minimum,  for  all  consistent  orderings. 

This  leaves  only  the  quantity  a^ ,  which  is  the  magnitude  of  the  pro- 

-(0) 

jection  of  the  initial  erroi  vector  £  onto  the  dominant  eigenvector, 
to  account  for  the  observed  behavior. 
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From  an  examination  of  equation  (20)  one  would  predict  that  if 
the  quantity  a^  were  in  fact  different  for  two  consistent  orderings 
then  the  error  vector  after  a  sufficient  number  of  iterations  would  be 
greater  or  less  for  the  one  scan  as  compared  to  the  other  by  a  fraction 
equal  to  the  ratio  of  the  two  a^'s.  Furthermore,  one  could  predict 
that  the  magnitude  of  the  error  vector  for  the  iteratively  slower  scan 
after  m  +  m1  iterations  would  be  equal  to  that  for  the  faster  scan 
after  m  iterations,  and  that  m'  would  be  independent  of  m.  Referring 
to  Fig.  10,  this  is  exactly  the  behavior  observed,  from  which  it  is 
concluded  that  the  quantity  aj,  corresponding  to  the  even-odd  scan 
is  less  than  the  aj  corresponding  to  the  standard  scan.  This  observa¬ 
tion  holds  for  all  the  problems  tested. 

That  the  aj’s  could  be  different  for  different  consistent  scans 
is  not  surprising  since,  although  all  of  the  eigenvalues  are  the  same, 
the  eigenvectors  are  known  to  be,  in  general,  different  for  different 
scans  (Ref  9:403),  and  hence  the  projection  of  the  initial  error  vector 
onto  the  dominant  eigenvector  will  in  general  be  differ3nt. 

The  data  indicate  that  the  even-odd  scan  is  preferable,  in  gener¬ 
al,  over  the  standard  scan  for  point  successive  overrelaxation.  The 
time  savings  achieved  will  vary  with  the  size  of  the  problem  and/or 
the  iteration  limit  required,  and  will,  be  reduced  asymptotically  to 
zero  as  the  number  of  iterations  increases  toward  infinity.  The  max¬ 
imum  savings  observed  for  few  iterations  was  on  the  order  of  15  to 
20$.  This  value  is  not  unreasonable  for  the  savings  to  bo  expected 
for  a  small  problem  (40  x  40  m^h)  with  a  required  error  vector 
reduction  of  three  or  four  orders  of  magnitude. 
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Cudahy's  Method 

One  of  the  initial  objectives  of  this  study  was  to  explain  why 
Cudahy's  method  was  faster  than  the  standard  scan  and  then  to  apply 
it,  or  a  variation  of  the  principal  with  block  overrelaxation,  to 
three  dimensions.  Unfortunately,  it  was  discovered  that  Cudahy's 
method  was  not  in  fact  asymptotically  faster  and  hence  the  expansion 
to  three  dimensions  was  not  attempted. 

From  the  data  presented  in  chapter  III,  Cudahy's  method  can  be 
seen  to  be  best  suited  for  the  transient  heat  equation  at  small  A  t. 
The  iterative  speed  of  his  method  can  be  seen  from  Tables  3.  and  5 
to  be  quite  favorable  when  compared  to  the  standard  scan  and  less  so 
when  compared  to  the  even-odd  scan.  The  data  in  Table  4,  for  A  t  = 
0.1,  include  the  results  for  Cudahy's  method  with  three  different 
relaxation  factors.  It  can  be  seen  for  an  iteration  limit  of  10“* 
that  the  highest  relaxation  factor  gives  the  greatest  iterative  speed. 
For  finer  convergence  criteria,  however,  the  CO  corresponding  to  the 
best  speed  decreases;  and  hence,  the  iterative  speed  is  a  function  of 
the  iteration  limit  as  well  as  the  relaxation  factor. 

The  method  developed  by  Cudahy  was  designed  to  accelerate  the 
convergence  of  point  successive  overrelaxation  by  scanning  the  bound¬ 
ary  values  into  the  mesh  as  quickly  possible.  Each  of  the  four 
scans  employed  takes  the  values  from  one  of  the  sides  of  the  mesh  and 
distributes  them  throughout  the  entire  mesh.  It  can  be  seen  from  Fig. 
8  that  this  produces  an  initial  rapid  convergence,  as  seems  reason¬ 
able,  but  that  the  method  is  eventually  iteratively  slower  than  the 
consistent  scans,  due  to  its  slower  asymptotic  rate  of  convergence. 
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Block  Overrel axatlon 

Although  this  present  work  was  confined  to  point  successive 
overrelaxation,  it  is  known  that  the  method  of  block  successive  over¬ 
relaxation  is  asymptotically  faster  by  approximately  40  per  oent 
(Ref  6:205). 


Much  of  the  theory  presented  in  chapter  II  for  point  successive 
overrelaxation  is  also  applicable  to  block  overrelaxation.  In  par¬ 
ticular  the  definition  Cf  and  techniques  for  finding  consistent 
orderings  are  identical.  It  is  only  necessary  to  partition  the  matrix 
A  of  equation  (l)  in  the  form 


A  = 


4l,l  ^1,2 

A  .  A . 

-2,1  -2,2 


A  4  A  _ 
-n,l  — n,2 


4,n 

A 

~2,n 


A 

-n,n 


(60) 


where  the  diagonal  submatrices  1  -  i  -  N,  are  square  and  non¬ 

singular.  The  associated  matrix  B  is  then  the  block  Jacobi  matrix  B 
defined  by  B  =  -D“^  A  +  I. 


As  a  practical  matter,  the  diagonal  submatrioes  of  A  of  equation 
(60)  must  be  tridiagonal  so  that  they  may  be  solved  directly. 

The  relevance  of  the  present  work  to  block  successive  overrelax- 
atlon  is  that  it  indicates  that  it  may  be  possible  to  find  a  consistent 
ordering  which  is  iteratively  faster  than  the  standard  method  of  line 
overrelaxation.  The  theory  does  not  specifically  exclude  this  possibil¬ 
ity;  however,  the  requirement  that  the  submatrices  be  tridiagonal  may 
make  the  search  difficult. 
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Recommendations 

Since  the  asymptotic  rate  of  convergence  for  point  and  block 
successive  overrelaxation  has  been  shown  by  Varga  to  be  optimum  for 
consistent  orderings  (Ref  6:125),  and  since  this  study  has  indicated 
that  the  iterative  speed  difference  between  two  consistent  orderings 
can  be  significant  for  small  problems,  the  author  feels  that  further 
research  into  accelerating  the  rate  of  convergence  could  be  profitably 
directed  toward  finding  that  consistent  ordering  which  yields  the 
greatest  initial  average  rate  of  convergence. 

Extension  of  the  method  to  three  dimensions  and  testing  of 
various  consistent  scans  in  that  domain  should  prove  worth-while. 

Other  iterative  methods,  such  as  the  alternating-direction 
implicit  methods  and  semi-iterative  methods,  have  been  developed  which 
appear  to  be  faster  than  point  successive  overrelaxation,  and  in  some 
cases,  block  successive  overrelaxation.  In  particular,  Varga  has 
reported  that  for  the  model  problem  (Laplace's  equation  in  a  rectangle) 
the  alternating-direction  implicit  methods  have  produced  convergence 
rates  as  much  as  38  times  as  great  as  for  point  successive  overrelax¬ 
ation  (Ref  6:228).  Varga  shows  that  these  methods  derive  their  speed 
primarily  from  the  use  of  different  accelerating  parameters  on  each 
iteration  (Ref  6:217).  A  significant  disadvantage  of  these  methods 
however,  in  addition  to  their  complexity,  is  the  fact  that  they  are 
not  rigorously  applicable  to  irregular  geometries. 

An  attempt  \ia.s  made  during  the  present  study  to  increase  the  rate 
of  convergence  for  the  even-odd  scam  by  the  use  of  non-stationary 
iterations  by  assigning  a  slightly  larger  than  optimum  relaxation 
factor  to  the  even  points  of  the  3can  and  a  slightly  smaller  tliar. 
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optimum  factor  for  the  odd  points  and  by  decreasing  the  difference 
for  subsequent  iterations,  but  the  method  proved  to  be  initially 
divergent  and  was  not  investigated  further.  It  may  be  possible  to 
develop  a  more  sophisticated  method  whioh  would  succeed. 

Use  of  the  techniques  employed  by  Cudahy,  combined  with  the  use 
of  different  relaxation  factors  on  each  iteration  might  lead  to 
greater  rates  of  convergence,  but  this  procedure  is  most  likely  to 
succeed  as  a  variation  of  the  alternating-direction  implicit  methods. 
Since  Cudahy's  method  has  the  ability  to  quickly  eliminate  large 
gradients  in  the  approximation  to  Laplace's  equation,  it  might  favor¬ 
ably  be  combined  in  an  iterative  procedure  with  the  even-odd  scan, 
wherein  Cudahy's  method  is  applied  initially  to  smooth  out  large 
gradients,  and  then  the  even-odd  scan,  with  its  greater  asymptotic 
rate  of  convergence,  is  used  to  obtain  the  desired  accuracy  through¬ 


out  the  mesh 
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